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In meta-analysis, between-study heterogeneity indicates the 
presence of effect-modifiers and has imphcations for the 
interpretation of resuhs in cost-effectiveness analysis and 
decision making. A distinction is usually made between 
true variability in treatment effects due to variation in 
patient populations or settings and biases related to the 
way in which trials were conducted. Variability in relative 
treatment effects threatens the external validity of trial evi- 
dence and limits the ability to generalize from the results; 
imperfections in trial conduct represent threats to internal 
validity. We provide guidance on methods for meta-regres- 
sion and bias-adjustment, in pairwise and network meta- 
analysis (including indirect comparisons), using illustrative 
examples. We argue that the predictive distribution of 
a treatment effect in a "new" trial may, in many cases, be 



more relevant to decision making than the distribution of 
the mean effect. Investigators should consider the relative 
contribution of true variability and random variation due 
to biases when considering their response to heterogeneity. 
In network meta-analyses, various types of meta-regression 
models are possible when trial-level effect-modifying cova- 
riates are present or suspected. We argue that a model with 
a single interaction term is the one most likely to be useful 
in a decision-making context. Illustrative examples of 
Bayesian meta-regression against a continuous covariate 
and meta-regression against "baseline" risk are provided. 
Annotated WinBUGS code is set out in an appendix. 
Key words: cost-effectiveness analysis; Bayesian meta- 
analysis; comparative effectiveness; systematic reviews. 
(Med Decis Making 2013;33:618-640) 



When combining results of different studies in 
a meta-analysis, heterogeneity (i.e., between- 
trials variation in relative treatment effects) has im- 
plications for the interpretation of results and deci- 
sion making. We provide guidance on techniques 
that can be used to explore the reasons for heteroge- 
neity, as recommended in the National Institute for 
Health and Clinical Excellence (NICE) Guide to 
Methods of Teclinology Appraisal.^ We focus 
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particularly on the implications of different forms 
of heterogeneity, the technical specification of mod- 
els that can estimate or adjust for potential causes of 
heterogeneity, and the interpretation of such models 
in a decision-making context. There is a considerable 
literature on the origins and measurements of het- 
erogeneity,^ and this is beyond the scope of this 
article. 

Heterogeneity in treatment effects is an indication 
of the presence of effect-modifying mechanisms — in 
other words, of interactions — ^between the treatment 
effect and the trial or trial-level variable. A distinc- 
tion is usually made between 2 kinds of interaction 
effects: clinical variation in treatment effects, result- 
ing from variation between treatment effects due to 
different patient populations, settings, or protocols 
across trials and deficiencies in the way the trial 
was conducted. The first type of interaction is said 
to represent a threat to the external validity of trials, 
limiting the extent to which one can generalize trial 
results from one situation to another. The second 
threatens the internal validity of trials: the trial deliv- 
ers a biased estimate of the treatment effect in its tar- 
get population, which may or may not be the same as 
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Figure 1 Posterior (solid) and predictive (dashed) densities for 
a treatment effect with mean = 0.7, standard deviation = 0.2, and 
heterogeneity (standard deviation) = 0.68. The area under the 
curve to the left of the vertical dotted line is the probability of a neg- 
ative value for the treatment effect. 

the target population for decision. Typically, these 
biases are considered to vary randomly over trials 
and do not necessarily have a zero mean. For exam- 
ple, the biases associated with markers of poor trial 
quality, such as lack of allocation concealment or 
lack of double blinding, have been shown to be asso- 
ciated with larger treatment effects.'^'* A general 
model for heterogeneity that encompasses both types 
of interaction has been proposed,^ but it is seldom 
possible to determine what the causes of heterogene- 
ity are and how much is due to true variation in clin- 
ical factors and how much is due to other unknown 
causes of biases. 

In a decision-making context, the response to high 
levels of heterogeneity is a critical issue. Investigators 
should compare the size of the treatment effect with 
the extent of between-trials variation. If the latter 
approximates the former, interpretation of mean 
treatment effects is difficult. ^'^ Figure 1 portrays a sit- 
uation in which a random effects [RE] model has been 
fitted. The posterior mean of the mean treatment 
effect is 0.70 with a posterior standard deviation 
(SD) of 0.2, making the mean effect clearly different 
from zero with a 95% credible interval (CrI) of 
(0.31, 1.09). However, the posterior mean of the 
between-trials standard deviation is ct = 0.68, compa- 
rable in size to the mean effect. What is a reasonable 



confidence interval for our prediction of the outcome 
of a future trial of infinite size? An approximate 
answer in classical statistics is found by adding the 
variance of the mean to the between-trials variance, 
which gives SD^ -I- ct^ = 0.50, giving a predictive stan- 
dard deviation of 0.71. Note that the 95% predictive 
interval is now (-0.69, 2.09), easily spanning zero 
effect, including a range of harmful effects. If we 
interpret these distributions in a Bayesian way, we 
would find that the probability that the mean effect 
is less than zero is only 0.0002, whereas the probabil- 
ity that a new trial would show a negative effect is 
much higher: 0.162 (Figure 1). 

This issue has been discussed before, and it has 
been proposed that, in the presence of heterogeneity, 
the predictive distribution, rather than the distribu- 
tion of the mean treatment effect, better represents 
our uncertainty about the comparative effectiveness 
of treatments in a future "rollout" of a particular 
intervention. In a Bayesian Markov chain Monte 
Carlo (MCMC) setting, a predictive distribution is 
easily obtained by drawing further samples from the 
distribution of effects: 

where d is the estimated (common) mean treatment 
effect and ct^ the estimated between-trial heterogene- 
ity variance. 

Where there are high levels of unexplained hetero- 
geneity, the implications of this recommendation for 
the uncertainty in a decision can be quite profound, 
and it is therefore important that the degree of hetero- 
geneity is not exaggerated and that the causes of the 
heterogeneity are investigated. The variance 
term in the predictive distribution should consist 
only of true variation between trial populations,^ 
but at present, there is no clear methodology to distin- 
guish between different sources of variation in treat- 
ment effect. Recent meta-epidemiological work on 
the determinants of between-study variation is begin- 
ning to shed some light on this.^^ 

Meta-regression is used to relate the size of a treat- 
ment effect obtained from a meta-analysis to potential 
effect modifiers — covariates — that may be character- 
istics of the trial participants or connected with the 
trial setting or conduct. These covariates may be cat- 
egorical or continuous. In common with other forms 
of meta-analysis, meta-regression can be based on 
aggregate (trial-level) outcomes and covariates, or 
individual patient data (IPD) and covariates when 
available. However, even if we restrict attention to 
randomized controlled trial (RCT) data, the study of 
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effect-modifiers is inherently observational^'^^ as it is 
not possible to randomize patients to one covariate 
value or another. As a consequence, meta-regression 
inherits all the difficulties of interpretation and infer- 
ence that attach to nonrandomized studies: con- 
founding, correlation between covariates, and, most 
important, the inability to infer causality from 
association. 

We provide guidance on methods for meta-regres- 
sion and bias-adjustment that can address the pres- 
ence of heterogeneity using aggregate data and trial- 
level covariates and discuss the limitations of this 
approach and the potential advantages of meta- 
regression with IPD. 



SUBGROUP EFFECTS (APPENDIX: EXAMPLE 1) 

We start with the simplest form of meta-regression: 
adjusting for subgroup effects in a pairwise meta- 
analysis. The theory set out in this section is then gen- 
eralized to other types of meta-regression and to mul- 
tiple treatments. The methods outlined here can be 
seen as extensions to the common generalized linear 
modeling framework for pairwise and network meta- 
analysis.^'^ 

In the context of treatment effects in RCTs, a sub- 
group effect can be understood as a categorical trial- 
level covariate that interacts with the treatment. 
The hypothesis would be that the size of the treat- 
ment effect is different, for example, in male and 
female patients or that it depends on age group or pre- 
vious treatment. Separate analyses can be carried out 
for each group and the relative treatment effects com- 
pared. However, if the models have random treat- 
ment effects, having separate analyses means 
having different estimates of between-trial variation. 
As there are seldom enough data to estimate the 
between-trial variation, it may make more sense to 
assume that it is the same for all subgroups. Further- 
more, running separate analyses does not immedi- 
ately produce the test of interaction that is required 
to reject the null hypothesis of equal effects. It is 
therefore preferable to have a single integrated analy- 
sis with a single between-trial variation term and an 
interaction term, |3, introduced on the treatment 
effect, as follows: 

^ik = M-i + {hih + (1) 

where Xy is the trial-level covariate for trial i, which 
can represent a subgroup, a continuous covariate, or 
baseline risk; is the linear predictor^'^ (e.g., the 



log-odds or mean) in arm k of trial i; |Xj are the trial- 
specific baseline effects in a trial i, treated as unre- 
lated nuisance parameters; 8j,ijt are the trial-specific 
treatment effects of the treatment in arm k relative 
to the treatment in arm 1 in that trial, with k=l,2 and 

. f 1 if u is true , , 

" \ 0 otlierwise ^ ' 

We can rewrite equation (1) as 

and note that the treatment and covariate interac- 
tion effects (8 and p) act only in the experimental 
arm, not in the control. For an RE model, the trial-spe- 
cific log-odds ratios come from a common distribu- 
tion: ?ii,i2~N{di2,(J^)- For a fixed effect (FE) model, 
we replace equation (1) with 0^^ = IJi'i + (c^i2 + 
(3xi)/{,5,^i}. In the Bayesian framework c/12. P and ct 
will be given independent (noninformative) priors: 
for example, c/12. P ~ N(0, 100^) and ct ~ Uniform(0, 5). 



META-REGRESSION MODELS IN NETWORK 
META-ANALYSIS 

In a network meta-analysis (NMA] context, vari- 
ability in relative treatment effects can also induce 
inconsistency^'* across pairwise comparisons. The 
methods introduced here are therefore also appropri- 
ate for dealing with inconsistency. Unless otherwise 
stated, when we refer to heterogeneity in this context, 
this can be interpreted as heterogeneity and/or 
inconsistency. 

A large number of meta-regression models can be 
proposed for NMA, each with very different implica- 
tions. There are 3 main approaches: unrelated inter- 
action terms for each treatment, exchangeable and 
related interaction terms, and one single interaction 
effect for all treatments. We argue that the third model 
is the most likely to have a useful interpretation in 
a decision-making context. 

Consider a binary between-trial covariate (e.g., pri- 
mary V. secondary prevention trials) in a case where s 
treatments are being compared. ^^'^^ We take treat- 
ment 1 (e.g., a placebo or standard treatment) as the 
reference treatment in the NMA. Following the 
approach to consistency models adopted previ- 
ously,^'^ we have (s - 1) basic parameters for the rela- 
tive treatment effects c/12, c/13, . . ., c/is of each 
treatment relative to treatment 1. The remaining (s — 
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l](s - 2)/2 treatment contrasts are expressed in terms 
of these parameters using the consistency equations: 
for example, the effect of treatment 4 compared with 
treatment 3 is written as c/34 = c?i4 — c/is.^^ We set out 
a range of fixed treatment effect interaction models, 
which are later extended to random treatment effects. 

1. Unrelated Treatment-Specific Interactions 

If there is an interaction effect between, say, pri- 
mary/secondary prevention and treatment, but these 
interactions are different for every treatment, the 
model will have as many interaction terms as there 
are basic treatment effects (e.g., ^12. P13, • • •, PiJ, 
each representing the additional (interaction) treat- 
ment effect in secondary prevention (compared 
with primary) in comparisons of treatments 2, 3, . . 
., s to treatment 1. These terms are exactly parallel 
to the main effects di2, di3, . . ., d-^g, which now repre- 
sent the treatment effects in primary prevention pop- 
ulations. As with the main effects, for trials 
comparing, say, treatments 3 and 4, the interaction 
term is the difference between the interaction terms 
on the effects relative to treatment 1, so that (334 = 
Pi4 - P13. The fixed treatment effects model for the 
linear predictor is 

= + {d\t,t - dua + (Piiij - (3) 

with tik representing the treatment in arm k of trial i, 
Xi the covariate/subgroup indicator, and I defined in 
equation (2). In all models, we set = Pn = 0. The 
remaining interaction terms are given unrelated 
vague prior distributions in a Bayesian analysis. 

The relative treatment effects in secondary preven- 
tion are c/12 + c?i3 + Pis, • ■ c/is + Pis- The 
interpretation is that the relative efficacy of each of 
the s treatments in primary prevention populations 
is entirely unrelated to their relative efficacy in sec- 
ondary prevention populations. 

2. Exchangeable and Related Treatment-Specific 
Interactions 

This model has the same structure as the model 
above, but now the (s — 1) "basic" interaction terms 
are drawn from a random distribution with a common 
mean and between-treatment variance: ^ij^~N{b,j^), 
for treatment k = 2, . . ., s. The mean interaction 
effect and its variance are estimated from the data, 
although informative priors that limit how similar 
or different the interaction terms are can be used. 



The interpretation is that there are real differences 
between the relative treatment effects, although these 
are centered on a common mean, with a common 
variance. 

3. Same Interaction Effect for All Treatments 

In this case, there is a single interaction term b that 
applies to the relative effects of all the treatments rel- 
ative to treatment 1 , so that p^^ = 6 for ic = 2 , . . . , s. The 
treatment effects relative to treatment 1, di2, di3...d-is 
in primary prevention, are all higher or lower by 
the same amount, b in secondary prevention: 
c/12 + b, di3 + b . . . dis + b. However, the effects 
of treatments 2, 3, . . ., s relative to each other in 
primary and secondary prevention populations are 
the same because the interaction terms cancel out. 
For example, consider the effect of treatment 4 rela- 
tive to treatment 3 in secondary prevention: 
du + b — {di3 + b) = di4 — c?i3, which is the same as 
in primary prevention. This means that the choice 
of reference treatment becomes important, and 
results for models with covariates are sensitive to 
this choice. 

Note that although we have presented the models 
in the context of subgroup effects, models for meta- 
regression with continuous covariates, including 
baseline risk, have the same structure. 

If an unrelated interactions model is being consid- 
ered (model 1), this requires 2 connected networks 
(one for each subgroup), including all the treatments, 
that is, with at least (s - 1) trials in each. With related 
interaction effects (model 2), it may not be necessary 
to have 2 complete networks, but a substantial num- 
ber of trials are needed to identify both random treat- 
ment effects and random interactions. However, to 
use such a model, there needs to be a clear rationale 
for exchangeability of interactions with a common 
mean and variance. One rationale could be to allow 
for different covariate effects for different treatments 
within the same class. Thus, treatment 1 is a standard 
or placebo treatment, whereas some of the treatments 
2, . . ., s belong to a "class." For example, one might 
imagine one set of exchangeable interaction terms 
for aspirin-based treatments for atrial fibrillation rel- 
ative to placebo and a second set of interactions for 
warfarin-based treatments relative to placebo. 

Although exchangeable and related interactions 
seem an attractive model in theory, they could 
lead to situations where the relative efficacy or the 
relative cost-effectiveness of a set of active treatments 
in the same class, and hence optimal treatment deci- 
sions, will depend on covariate values. However, 
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differences between interaction terms are unlikely to 
be robust, and treatment decisions could then be 
driven by statistically insignificant interaction terms. 
Therefore, although these more complex models can 
be fitted, their best role may be in exploratory 
analyses. Here we explore only model 3, which 
assumes an identical interaction effect across all 
treatments with respect to the reference treatment. 

Network Meta-Regression with a Continuous 
Covariate 

When dealing with a continuous covariate, we will 
use centered covariate values. This is achieved by 
subtracting the mean covariate value, x, from each x,-. 
To fit a random treatment effects model that assumes 
a common interaction effect for all treatments with 
covariate centering, we rewrite equation (3) as 

= ^^i + {Kik + - (4) 
where Pn = 0, = b [k = 2, . . . ,7] and bijk~N 

We retain the treatment-specific interaction effects 
but set them all equal to b, so that the terms cancel out 
in comparisons not involving the reference treat- 
ment. To produce treatment effect estimates at cova- 
riate value z, the outputs from this model should be 
"uncentered": dik — b{x — z), k = 2, . . ., s. 

Example: Certolizumab Meta-Regression on Mean 
Disease Duration (Appendix: Example 2) 

Data are available from a review of trials of certolizu- 
mab pegol (CZP) for the treatment of rheimiatoid arthri- 
tis in patients who had failed on disease-modifying 
antirheumatic drugs, including methotrexate (MTX).^° 
Twelve MTX controlled trials were identified, compar- 
ing 6 different treatments with placebo (Figure 2). 

Table 1 shows the number of patients who have 
improved by at least 50% on the American College 
of Rheumatology scale (ACR50) at 6 months, r^^, out 
of all included patients, rtu^, for each arm of the 
included trials, along with the mean disease duration 
in years for patients in each trial, Xj (i = 1, . . ., 12; ic= 1, 
2). The crude odds ratios (ORs) from Table 1 are plot- 
ted (on a log-scale] against mean disease duration in 
Figure 3, with the numbers 2 to 7 representing the 
OR of that treatment relative to placebo plus MTX 
(chosen as the reference treatment). Example 1 in 
the Appendix gives details of priors and the Win- 
BUGS^^ code for the model in equation (4). 

Table 2 shows the results of fitting fixed and ran- 
dom treatment effects NMA^'^ and interaction models 



CZP + MTX (2) 






Adalimumab+ MTX (3) 


/ }^ 
Placebo + MTX (1) ^CT''''^ 


^ — Etanercept+ MTX (4) 


\ 

\2 ^\ 


Infliximab+ MTX (5) 


\ \ Rituximab+ MTX (6) 
Tocilizumab+ MTX (7) 


Figure 2 Certolizumab example^": 
necting 2 treatments indicate that 


treatment network. Lines con- 
a comparison between these 



treatments has been made. The numbers on the lines indicate 
how many randomized controlled trials compare the 2 connected 
treatments. CZP, certolizumab pegol; MTX, methotrexate. 

with disease duration as the covariate. The estimated 
ORs for different durations of disease are represented 
by the parallel lines in Figure 3. Note that the treat- 
ment effects obtained are the estimated log-odds 
ratios at the mean covariate value [x = 8.21 years in 
this case) and that credible interval lines have not 
been included in Figure 3. 

The deviance information criterion (DIC) provides 
a measure of model fit that penalizes model complex- 
ity^^ — lower values of the DIC suggest a more parsi- 
monious model. We calculate the DIC as the sum of 
the posterior mean of the residual deviance, Dres, 
and the leverage, po-^^ The DIC and posterior means 
of the residual deviances in Table 2 do not decisively 
favor a single model (differences less than 3 or 5 are 
not considered meaningful). The fit of the FE model 
is improved by including the covariate interaction 
term b, which also has a CrI that does not include 
zero. Addition of the covariate to the RE model 
reduces the heterogeneity from a posterior median 
of 0.34 to 0.28, but the CrI for the interaction parame- 
ter b includes zero. In this example, the meta-regres- 
sion models are all reasonable but not strongly 
supported by the evidence, although the finding of 
smaller treatment effects with a shorter disease dura- 
tion has been reported previously. 

Should the use of biologies be confined to patients 
whose disease duration was above a certain thresh- 
old? This seems reasonable, but it would be difficult 
to determine this threshold on the basis of this evi- 
dence alone. The slope is largely determined by treat- 
ments 3 and 7 (adalimumab and tocilizumab), which 
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Table 1 Certolizumab Example^": Number of Patients Achieving ACR50 at 6 Months, r, Out of the Total 
Number of Patients, n, in Arms 1 and 2 of the 12 Trials and Mean Disease Duration (in Years) for Patients in 

Trial i, Xi 



Arm 1 Arm 2 



Study Name 


Treatment in 
Arm 1, fii 


Treatment in 
Arm 2, (,2 


No. Achieving 
ACR50, Til 


Total No. of 
Patients, iiji 


No. Achieving 
ACR50, 


Total No. of 
Patients, n;2 


Mean Disease 
Duration (Years), X; 


RAPID 1 


Placebo 


CZP 


15 


199 


146 


393 


6.15 


RAPID 2 


Placebo 


CZP 


4 


127 


80 


246 


5.85 


Kim 2007 


Placebo 


Adalimumab 


9 


63 


28 


65 


6.85 


DE019 


Placebo 


Adalimumab 


19 


200 


81 


207 


10.95 


ARMADA 


Placebo 


Adalimumab 


5 


62 


37 


67 


11.65 


Weinblatt 1999 


Placebo 


Etanercept 


1 


30 


23 


59 


13 


START 


Placebo 


Infliximab 


33 


363 


110 


360 


8.1 


ATTEST 


Placebo 


Infliximab 


22 


110 


61 


165 


7.85 


Abe 2006*' 


Placebo 


Infliximab 


0 


47 


15 


49 


8.3 


Strand 2006 


Placebo 


Rituximab 


5 


40 


5 


40 


11.25 


charisma'^ 


Placebo 


Tocilizumab 


14 


49 


26 


50 


0.915 


OPTION 


Placebo 


Tocilizumab 


22 


204 


90 


205 


7.65 



All trial arms had methotrexate in addition to the placebo or active treatment. CZP, certolizumab pegol. 
'^Three months was used when this was not available. 
•"ACRSO at 3 months. 



are the only treatments trialed at more than 1 disease 
duration and appear to have different effects at each 
duration (Figure 3). However, the linearity of rela- 
tionships is highly questionable, and the prediction 
of negative effects for treatment 6 (rituximab) is not 
plausible. This suggests that other explorations of 
the causes of heterogeneity in this network should 
be undertaken (see below). 

Network Meta-Regression on Baseline Risk 

The meta-regression model presented above can be 
extended to use the baseline risk in each trial as 
a covariate, taking into account the error in the esti- 
mation of baseline risk and its correlation to the OR. 
The model is the same as in equation (4) but now Xi 
= |Xy, the trial-specific baseline for the control arm in 
each trial. An important property of this Bayesian for- 
mulation is that it takes the "true" baseline (as esti- 
mated by the model) as the covariate and 
automatically takes the uncertainty in each into 
account. Naive approaches, which regress 
against the observed baseline risk, fail to take the cor- 
relation between the treatment effect and baseline 
risk, as well as the consequent regression to the 
mean phenomenon, into account. 

Example: Certolizumab Meta-Regression on Base- 
line Risk (Appendix: Example 3) 



Both FE and RE models with a common interaction 
term were fitted. Due to covariate centering, the treat- 
ment effects in Table 3 for the models with covariate 
adjustment are the effects for patients with a baseline 
log-odds ACR50 of —2.421, the mean of the observed 
log-odds on treatment 1. This corresponds to a base- 
line ACR50 probability of 0.082. 

In both models, the 95% Cris for covariate term 
exclude zero, suggesting a strong interaction effect 
between the baseline risk and the treatment effects. 
The estimated RE model is shown in Figure 4: the dif- 
ferences between the lines represent treatment differ- 
ences controlling for baseline risk. However, note 
that there is considerable uncertainty around these 
lines, which is not represented in Figure 4. The DIG 
statistics and the posterior means of the residual devi- 
ance also marginally favor the RE model with the 
covariate. In interpreting the results of this model, it 
should be noted that treatment effect is being defined 
on a percent change scale: the extent of change is 
therefore automatically increased as the baseline is 
lowered. Hence, the observed interaction may to 
some extent be an artifact of the measurement scale. 



BIAS AND BIAS-ADJUSTMENT 

If it is thought that heterogeneity is due to bias 
in some of the included studies, models for bias- 
adjustment can be used. The difference between 
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Figure 3 Certolizumab example^" : plot of the crude odds ratio ( OR) 
(on a log-scale) of the 6 active treatments relative to placebo plus 
methotrexate (MTX) against mean disease duration (in years). 
For plotting purposes, the odds of response on placebo plus MTX 
for the Abe 2006 study were assumed to be 0.01. The plotted numbers 
refer to the treatment being compared with placebo plus MTX, and 
the lines represent the relative effects of the following treatments 
compared with placebo plus MTX based on a random effects meta- 
regression model (from top to bottom): etanercept plus MTX (treat- 
ment4, dotted green line) , certolizumab pegol (CZP) plus MTX (treat- 
ment 2, solid black line), tocilizumab plus MTX (treatment 7, short- 
longdash purple line), adalimumab plus MTX (treatment 3, dashed 
red line), infliximab plus MTX (treatment 5, dot-dashed dark blue 
line), and rituximab plus MTX (treatment 6, long-dashed black 
line). Odds ratios above 1 favor the plotted treatment, and the hori- 
zontal line (thin dashed) represents no treatment effect. 

"bias-adjustment" and the meta-regression models 
described above is slight but important. In meta-regres- 
sion, we concede that even within the formal scope of 
the decision problem, there are distinct differences in 
relative treatment efficacy. In bias-adjustment, we 
have in mind a target population for decision making, 
but the evidence available, or at least some of the evi- 
dence, provides biased, or potentially biased, estimates 
of the target parameter, perhaps because the trials have 
internal biases, they concern different populations or 
settings, or both. 

Adjustment for Bias Based on Meta-Epidemiologi- 
cal Data 

Confronted by trial evidence of mixed quality, 
investigators have 3 options: 1] to restrict attention 



to studies of high quality, ignoring what may be a sub- 
stantial proportion of the evidence; 2) to include tri- 
als of both high and low quality in a single analysis, 
which risks delivering a biased estimate of the treat- 
ment effect; or 3) to use all the data but simulta- 
neously adjust and down-weight the evidence from 
studies at risk of bias. The latter uses information 
on the expected bias, as well as between-study vari- 
ation in bias, derived from statistical analysis of 
meta-epidemiological data. The analysis assumes 
that the study-specific biases in the data set of 
interest are exchangeable with those in the meta- 
epidemiological data used to provide the prior dis- 
tributions used for adjustment. This approach can 
form the basis for sensitivity analyses to show 
whether the presence of studies at risk of bias, with 
potentially over-optimistic results, has an impact 
on the decision. 

It is expected that this bias-adjustment method 
will be used more when more data are available on 
how the degree of bias depends on the nature of the 
outcome measure and the clinical condition.* In prin- 
ciple, the same form of bias-adjustment could be 
extended to other types of bias or to mixtures of 
RCTs and observational studies. Each of these exten- 
sions, however, depends on the applicability of anal- 
yses of very large meta-epidemiological data sets that 
are starting to become available. ^^'^^ 

Estimation and Adjustment for Bias in Network 
Meta-Analysis 

In NMA, if we assume that the mean and variance 
of the study-specific biases are the same for each 
treatment, it is possible to simultaneously estimate 
the treatment effects and the bias effects in a single 
analysis and to produce treatment effects that are 
based on the entire body of data and also adjusted 
for bias.^^ If A is placebo or standard treatment and 
B, C, and D are all active treatments, it would be rea- 
sonable to expect the same bias distribution to apply 
to the AB, AC, and AD trials. But it is less clear what 
the direction of the bias should be in BC, BD, and CD 
trials. Assuming that the average bias is always in 
favor of the newer treatment, this becomes a model 
for novelty bias.^^'^^ Another approach might be to 
propose a separate mean bias term for active v. active 
comparisons.^^ 

This method can in principle be extended to 
include syntheses that are mixtures of trials and 
observational studies, and it can also be extended to 
any form of "internal" bias. A particularly interesting 
application is to "small-study bias," which is one 
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Table 2 Certolizumab Example : Posterior Mean, SD, and 95% CrI for the Interaction Estimate [b) and Log- 
Odds Ratios dxY of Treatment Y relative to Treatment X 



No Covariate Covariate "Disease Duration" 



FE RE" EE RE" 





Mean 


SD 


CrI 


Mean/ Median 


SD 




CrI 




Mean 


SD 


CrI 


Mean/ Median 


SD CrI 


b 


NA 


NA 


NA 


NA 


NA 




NA 




0.14 


0.06 


(0.01, 0.26) 


0.14 


0.09 (-0.03, 0.32) 




2.21 


0.25 


(1.73, 2.72) 


2.27 


0.39 


(1. 


.53, 3. 


.10) 


2.50 


0.29 


(1.96, 3.08) 


2.57 


0.42 (1.79, 3.44) 


C?13 


1.93 


0.22 


(1.52, 2.37) 


1.97 


0.33 


(1. 


.33, 2. 


.64) 


1.66 


0.25 


(1.19, 2.16) 


1.71 


0.34 (1.04, 2.41) 




3.47 


1.34 


(1.45, 6.74) 


3.46 


1.41 


(1. 


.26, 6. 


.63) 


2.82 


1.34 


(0.71, 5.96) 


2.77 


1.42 (0.42, 6.01) 


di5 


1.38 


0.17 


(1.06, 1.72) 


1.48 


0.33 


(0. 


.90, 2. 


.21) 


1.40 


0.17 


(1.08, 1.74) 


1.48 


0.30 (0.95, 2.15) 


dw 


0.00 


0.71 


(-1.40, 1.39) 


0.01 


0.82 ( 


-1, 


.61, 1. 


.63) 


-0.42 


0.73 


(-1.86, 1.04) 


-0.44 


0.84 (-2.08, 1.21) 


di7 


1.65 


0.22 


(1.22, 2.10) 


1.56 


0.38 


(0. 


.77, 2. 


.28) 


1.98 


0.28 


(1.45, 2.53) 


2.00 


0.45 (1.12, 2.93) 


a 


NA 


NA 


NA 


0.34 


0.20 


(0. 


.03, 0. 


.77) 


NA 


NA 


NA 


0.28 


0.19 (0.02,0.73) 


D ^ 

■'-'res 


37.6 






30.9 










33.8 






30.2 




Pd 


18.0 






21.2 










19.0 






21.3 




DIG 


55.6 






52.1 










52.8 






51.4 





Posterior median, standard deviation (SD), and 95% credible interval (Crl) of the between-trial heterogeneity (cr) for the number of patients achieving 
ACR50 for the fixed effects (FE) and random effects (RE) models with and without covariate "disease duration" and measures of model fit: posterior 
mean of the residual deviance (-Dres), effective number of parameters [po), and deviance information criterion (DIG). Results are based on 100,000 iterations 
from 3 independent chains after a burn-in of 40,000 iterations. Treatment codes are given in Figure 2. NA, not applicable. 
^Using informative prior for a (for details, see Appendix). 
^Compare with 24 data points. 



Table 3 Certolizumab Example^^: Posterior Mean, SD, and 95% CrI for the Interaction Estimate [b) and Log- 
Odds Ratios dxY of Treatment Y Relative to Treatment X 





FE 






RE 




Mean 


SD 


CrI 


Mean/Median 


SD 


CrI 



b 


-0.93 


0.09 


(-1.03,-0.69) 


-0.95 


0.10 


(-1.10, -0.70) 


d-i_2 


1.85 


0.10 


(1.67, 2.06) 


1.83 


0.24 


(1.35, 2.29) 


di3 


2.13 


0.11 


(1.90, 2.35) 


2.18 


0.22 


(1.79, 2.63) 


di4 


2.08 


0.34 


(1.47, 2.80) 


2.04 


0.46 


(1.19, 2.94) 


di5 


1.68 


0.10 


(1.49, 1.86) 


1.71 


0.22 


(1.30, 2.16) 


die 


0.36 


0.50 


(-0.72, 1.27) 


0.37 


0.59 


(-0.86, 1.45) 


di7 


2.20 


0.14 


(1.93, 2.46) 


2.25 


0.27 


(1.75, 2.79) 


(J 


NA 


NA 


NA 


0.19 


0.19 


(0.01, 0.70) 


^res 


27.3 






24.2 






Pd 


19.0 






19.4 






DIG 


46.3 






43.6 







Posterior median, standard deviation [SD), and 95% credible interval (CrI) of the betv^^een-trial heterogeneity (ct) for the number of patients achieving 
ACR50 for the fixed effects (FE) and random effects (I^E) models with covariate "baseline risk" and measures of model fit: posterior mean of the residual 
deviance [Dres), number of parameters (p^). deviance information criterion (DIG). Results are based on 100,000 iterations from 3 independent chains 
after a burn-in of 60,000 iterations. Treatment codes are given in Figure 2. NA, not applicable. 
'^Compare with 24 data points. 



interpretation of "publication bias." The idea is that 
smaller studies have greater bias. The "true" treat- 
ment effect can therefore be conceived as the effect 
that would be obtained in a study of infinite size, 
which is taken to be the intercept in a regression of 
the treatment effect against the study variance. '^"'^^ 



Elicitation of Bias Distributions 

This method^^ is conceptually the simplest of all 
bias-adjustment methods, applicable to trials and 
observational studies alike, but it is also the most dif- 
ficult and time-consuming. One of its advantages is 
that it can be used when the number of trials is 
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baseline odds of response (log-scale) 



Figure 4 Certolizumab example^" : plot of the crude odds ratio ( OR] 
of the 6 active treatments relative to placebo plus methotrexate 
(MTX) against odds of baseline response on a log-scale. For plot- 
ting purposes, the odds of response on placebo plus MTX for the 
Abe 2006 study were assumed to be 0.01. The plotted numbers refer 
to the treatment being compared with placebo plus MTX, and the 
lines represent the relative effects of the following treatments 
(from top to bottom) compared with placebo plus MTX based on 
a random effects meta-regression model: tocilizumab plus MTX 
(treatment 7, short-long dash purple line], adalimumab plus MTX 
(treatment 3, dashed red line), etanercept plus MTX (treatment 4, 
dotted green], certolizumab pegol (CZP) plus MTX (treatment 2, 
solid black line], infliximab plus MTX (treatment 5, dot-dashed 
dark blue line), and rituximab plus MTX ( treatment 6, long-dashed 
black line]. Odds ratios above 1 favor the plotted treatment, and the 
horizontal line ( dashed ) represents no treatment effect. 

insufficient for meta-regression approaches. Each 
study is considered by several independent experts 
using a predetermined protocol that itemizes a series 
of potential internal and external biases. Each expert 
is asked to provide information that is used to develop 
a combined bias distribution. The mean and variance 
of the bias distributions are statistically combined 
with the original study estimate and its variance to cre- 
ate a new, adjusted estimate of the treatment effect in 
each study. The final stage is a conventional synthesis 
of the adjusted study-specific estimates using stan- 
dard pairwise or network meta-analysis methods. 

DISCUSSION 

We have detailed methods to explore or explain 
heterogeneity in pairwise and network meta-analyses 



using study-level covariates, including baseline risk 
(variation in "baseline" natural history is dealt with 
in another tutorial in this series'^'^), although we 
have not covered the closely related topic of outlier 
detection."'^*'^^ 

However, meta-regression based on aggregate data 
and study-level covariates suffers from several bias 
and confounding problems as well as typically hav- 
ing low power. When a categorical covariate is con- 
sidered, one can contrast a within-trial comparison 
(e.g., treatment effects reported separately for males 
and females) and a between-trial comparison where 
different trials are run on male and female patients. 
The contrast is similar to the difference between 
a paired and an unpaired t test. With between-trial 
comparisons, a given covariate effect (i.e., interac- 
tion) will be harder to detect as it has to be distin- 
guishable from the "random noise" created by the 
between-trial variation. However, for within-trial 
comparisons, the between-trial variation is con- 
trolled for, and the interaction effect needs only to 
be distinguishable from sampling error. A further 
drawback of between-trial comparisons is that, 
because the number of observations (trials) may be 
very low while the precision of each trial may be rel- 
atively high, it is quite possible to observe a highly 
statistically significant relation between the treat- 
ment effect and the covariate that is entirely spuri- 
ous. '^^ Furthermore, between-trial comparisons are 
more vulnerable to ecologic bias or ecologic fallacy, 
where, for example, a linear regression coefficient of 
treatment effect against the covariate in the 
between-trial case can be entirely different from the 
coefficient for the within-trial data. 

With continuous covariates and IPD, not only does 
the within-trial comparison avoid ecological bias, but 
it also has far greater statistical power to detect a true 
covariate effect. This is because the variation in 
patient covariate values will be many times greater 
than the variation between the trial means, ^'''^^ and 
the precision in any estimated regression coefficient 
depends directly on the variance in covariate values. 
For these reasons, IPD meta-analyses and meta- 
regression are taken to be the "gold standard, "''^ 
although IPD is not available in most situations. 

Two broad approaches to IPD meta-regression 
have been considered."*^ In the 2-step approach, the 
analyst first estimates the interaction effect sizes 
and their standard errors from each study and then 
conducts a standard meta-analysis on these summa- 
ries. This is appropriate for inference on the existence 
of an interaction, but for decision making, we recom- 
mend a 1-step method in which main effects and 
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interaction are estimated simultaneously. This is 
because the parameter estimates of the main effects 
and interaction terms will be correlated, and their 
joint uncertainty can most easily be propagated 
through the decision model by estimating them simul- 
taneously. IPD RE pairwise meta-analysis models have 
been developed for continuous binary** sur- 
vival,*^ and ordinal*^ variables. Although most of the 
models are presented in the single pairwise compari- 
son context, it is possible to extend them to an NMA 
context. *^'*^~^^ Criteria for determining the potential 
benefits of IPD to assess patient-level covariates have 
recently been outlined. 

When IPD is available only in a subset of studies, it 
is possible to incorporate both types of data in the 
same analysis using dual effect models. *^'^°'^^~^^ 
This makes the best possible use of all available 
data, but a decision has to be made on whether 
between-study variability is to be included in the esti- 
mation of effects. Statistical tests are unlikely to have 
sufficient power to inform a decision. In an NMA 
with IPD, we would recommend use of models with 
a single interaction term for each covariate, at least 
within a class of treatments, for decision making, 
but more complex structures have been attempted.*^ 

Finally, it needs to be appreciated that in cases 
where the covariate does not interact with the treat- 
ment effect but modifies the baseline risk, the effect 
of pooling data over the covariate is to bias the esti- 
mated treatment effect toward the null effect. This 
is a form of ecologic bias known as aggregation 
bias,^^ which does not affect strictly linear models, 
where pooling data across such covariates will not 
create bias. Usually, it is significant only when both 
the covariate effect on baseline risk and the treatment 
effect are quite strong. It is a particular danger in 



survival analysis because the effect of covariates 
such as age on cancer risk can be particularly marked 
and because the log-linear models routinely used are 
highly nonlinear. When covariates that affect risk are 
present, even if they do not modify the treatment 
effect, the analysis must be based on pooled estimates 
of treatment effects from a stratified analysis for 
group covariates and regression for continuous cova- 
riates and not on treatment effects estimated from 
pooled data. 

Bias-adjustment methods are another form of 
accounting for differences between studies. These 
methods should be considered semi-experimental 
as there is a need for further experience with applica- 
tions and for further meta-epidemiological data on 
the relationships between the many forms of internal 
bias that have been proposed. However, they repre- 
sent reasonable and valid methods for bias-adjust- 
ment and are likely to be superior to no bias- 
adjustment in situations where data are of mixed 
quality. At the same time, caution is required as the 
method is essentially a meta-regression based on 
"between-studies" comparisons, with no certainty 
of a causal mechanism. 

We have suggested that choice between models 
with and without the covariate or bias-adjustment 
coefficients should be based on examining the DIG 
and the CrI of the regression coefficient. However, 
in general, the DIG is not able to inform choice 
between RE models with and without covariates as 
both will tend to fit equally well and have a similar 
effective number of parameters. Other considera- 
tions, such as a reduction in the heterogeneity, will 
play a greater role in choosing between RE models. 
For FE models with and without covariates, the DIG 
is suitable for model choice. 



APPENDIX 

Heterogeneity and Meta-Regression — Illustrative Examples and WinBUGS Code 

This appendix gives illustrative WinBUGS code for all the examples presented in the main article. All program- 
ming code is fully annotated. The program codes are printed here but are also available as WinBUGS system files 
from http://www.nicedsu.org.uk. Users are advised to download the WinBUGS files from the website instead of 
copying and pasting from this document. We have provided the codes as complete programs. However, the major- 
ity of each meta-regression program is identical to the programs in Dias and others. ^^'^^ We have therefore high- 
lighted the main differences in bold, to emphasize the modular nature of the code. The code presented is 
completely general and will be suitable for fitting pairwise or network meta-analyses with any number of treat- 
ments and multiarm trials. We also provide an indication of the relevant parameters to monitor for inference and 
model checking for the various programs. The nodes to monitor for the FE models are the same as those for the 
RE models, except that there is no heterogeneity parameter. 
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Table Al gives an index of the programs and their relation to the examples in the main text. Note that for 
each example, there are RE and FE versions of the code. All FE code can be run using the same data structure 
described for RE. 



Table Al Index of WinBUGS Code with Details of Examples and Sections Where They Are Described 





Program 


Fixed or Random Effects 


Example Name 




Model Specification 


1 


(a) 


RE 


Statins (subgroups) 


Meta-reg 


jression with subgroups 




(b) 


FE 








2 


(a) 


RE 


Certolizumab (mean disease dura- 


Meta-reg 


jression with continuous covariate 




(b) 


FE 


tion) 






3 


(a) 


RE 


Certolizumab (basebne risk) 


Meta-reg 


;ression with adjustment for base- 




(b) 


FE 




line risk 





FE, fixed effects; RE, random effects. 



Example 1. Statins: Meta-Regression with Subgroups 

A meta-analysis of 19 trials of statins for cholesterol lowering v. placebo or usual care^^ included some 
trials on which the aim was primary prevention (patients included had no previous heart disease) and others 
on which the aim was secondary prevention (patients had previous heart disease). Note that the subgroup 
indicator is a trial-level covariate. The outcome of interest was all-cause mortality, and the data are presented 
in Table A2. The potential effect-modifier, primary v. secondary prevention, can be considered a subgroup in 
a pairwise meta-analysis of all the data, or 2 separate meta-analyses can be conducted on the 2 types of study. 

Table A2. Meta- Analysis of Statins v. Placebo for Cholesterol Lowering in Patients with and without 
Previous Heart Disease^^: Number of Deaths Due to All-Cause Mortality in the Control and Statin Arms of 

19 Randomized Controlled Trials 



Placebo/Usual Care Statin 



Trial ID 


No. of Deaths, r,j 


No. of Patients, n,j 


No. of Deaths, r^^ 


No. of Patients, n^^ 


Type of Prevention, 


1 


256 


2223 


182 


2221 


Secondary 


2 


4 


125 


1 


129 


Secondary 


3 


0 


52 


1 


94 


Secondary 


4 


2 


166 


2 


165 


Secondary 


5 


77 


3301 


80 


3304 


Primary 


6 


3 


1663 


33 


6582 


Primary 


7 


8 


459 


1 


460 


Secondary 


8 


3 


155 


3 


145 


Secondary 


9 


0 


42 


1 


83 


Secondary 


10 


4 


223 


3 


224 


Primary 


11 


633 


4520 


498 


4512 


Secondary 


12 


1 


124 


2 


123 


Secondary 


13 


11 


188 


4 


193 


Secondary 


14 


5 


78 


4 


79 


Secondary 


15 


6 


202 


4 


206 


Secondary 


16 


3 


532 


0 


530 


Primary 


17 


4 


178 


2 


187 


Secondary 


18 


1 


201 


3 


203 


Secondary 


19 


135 


3293 


106 


3305 


Primary 



The number of deaths in arm k of trial i, r^, is assumed to have a binomial hkelihood r.^ ~ Binomial(/7;^^,«;^), 
1 = 1,..., 19; k= 1,2. Defining x- as the trial-level subgroup indicator such that 




0 if study i is a primary prevention study 

1 if study i is a secondary prevention study 
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our interaction model is 0 = \i. + . j^, + [3x. j , where 0 =logit( j!?.^ ) is the linear predictor (see Dias 
and others^^). In this setup, |i- represent the log-odcis of the outcome in the "control" treatment (i.e., the treat- 
ment indexed 1), and 5, j2 are the trial-specific log-odds ratios of success on the treatment group compared 
with control for primary prevention studies. See the main text for further details of the model. 

WinBUGS code to fit 2 separate fixed or random effects models is given in the appendix to Dias and oth- 
ers." The WinBUGS code for a single analysis with an interaction term for subgroup with RE is given in 
program 1(a), and the FE code is given in program 1(b). Although this example includes only 2 treatments, 
the code presented below can also be used for subgroup analysis with multiple treatments and including 
multiarm trials. 



Program 1(a): 
Example) 



Rinomial Likelihood, Logit Link, Random Effects, Meta-Regression with Subgroups (Statins 



# Binomial likelihood, logit link, subgroup 

# Random effects model for multi-arm trials 
modelj 

for(i in l:ns]( 
w[i,l] <- 0 
delta[i,l] <- 0 
mu[i] ~ dnorm(0,.0001) 
for (k in l:na[i]) | 



#*** PROGRAM STARTS 

# LOOP THROUGH STUDIES 

# adjustment for multi-arm trials is zero for control arm 

# treatment effect is zero for control arm 

# vague priors for all trial baselines 

# LOOP THROUGH ARMS 

# binomial likelihood 



r[i,k] ~ dbin(p[i,k],n[i,k]) 
# model for linear predictor, covariate effect relative to treat in arm 1 
logit(p[i,k]] <- mu[i] -i- deha[i,k] -i- (beta[t[i,k]]-beta[t[i,l]]) * x[i] 

rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]] - log(n[i,k]-rhat[i,k]))) 



resdev[i] <- sum(dev[i,l:na[i]]) 

for (k in 2:na[i]) | 

delta[i,k] ~ dnorm(md[i,k],taud[i,k]) 
md[i,k] <- d[t[i,k]] - d[t[i,l]] -i- sw[i,k] 
taud[i,k] <- tau *2*(k-l)/k 
w[i,k] <- (delta[i,k] - d[t[i,k]] -I- d[t[i,l]]) 
swr[i,k] <- sum(w[i,l:k-l])/(k-l) 



# summed residual deviance contribution for this trial 

# LOOP THROUGH ARMS 

# trial-specific LOR distributions 

# mean of LOR distributions (with multi-arm trial correction) 

# precision of LOR distributions (with multi-arm trial correction) 

# adjustment for multi-arm RCTs 

# cumulative adjustment for multi-arm trials 



totresdev <- sum(resdev[]) 
d[l]<-0 
beta[l] <- 0 
for (k in 2:nt)| 

d[k] ~ dnorm(0,.0001) 

beta[k] <- B 

I 

B ~ dnorm(0,.0001) 

sd ~ dunif(0,5) 
tau <- pow(sd,-2) 
1 



# Total Residual Deviance 

# treatment effect is zero for reference treatment 

# covariate effect is zero for reference treatment 

# LOOP THROUGH TREATMENTS 

# vague priors for treatment effects 

# common covariate effect 

# vague prior for covariate effect 

# vague prior for between-trial SD 

# between-trial precision = (1 /between-trial variance) 

# *** PROGRAM ENDS 



To obtain posterior summaries for other parameters of interest, the nodes d, B, and sd need to be moni- 
tored. To obtain the posterior means of the parameters required to assess model fit and model comparison, 
dev, totresdev, and the DIG (from the WinBUGS DIG tool) need to be monitored. 
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Additional code can be added before the closing brace to estimate all the pairwise log-odds ratios and odds 
ratios and to produce estimates of absolute effects, given additional information on the absolute treatment effect 
on one of the treatments, for given covariate values. For further details on calculating other summaries from the 
results and on converting the summaries onto other scales, refer to the appendix in Dias and others. 

################################################################################ 

# Extra code for calculating all odds ratios and log odds ratios, and absolute effects, for covariate 

# values in vector z, with length nz (given as data) 

################################################################################ 
for (k in lmt){ 

for (j in l:nz) { dz[i,k] <- d[k] + (beta[k]-beta[l])*z[j] } # treatment effect when covariate = z[j] 

} 

# pairwise ORs and LORs for all possible pair-wise comparisons 
for (c in l:(nt-l)) { 

for (k in (cH-l):nt) { 

# when covariate is zero 

or[c,k] <- exp(d[k] - d[c]] 
lor[c,k] <- (d[k]-d[c]) 

# at covariate=z[j] 

for (j in l:nz) { 

orz[j,c,k] <- exp(dz[j,k] - dz[j,c]) 
lorz[j,c,k] <- (dz[j,k]-dz[j,c]] 

} 

} 

} 

# Provide estimates of treatment effects T[k] on the natural (probability) scale 

# Given a Mean Effect, meanA, for 'standard' treatment 1, with precision (1/variance) precA, and covariate 
value z[]] 

A ~ dnorm(meanA,precA) 
for (k in Imt) { 
for (j in l:nz){ 

logit(T[j,k]) <- A -H d[k] + (beta[k]-beta[l]) * z[j] 

} 

} 

For a meta-regression with 2 subgroups, vector z would be added to the list data statement as list(z=c(l), 
nz=l). 

The data structure is identical to that presented in Dias and others^^ but now has an added column x[] that 
represents the value of the covariate (taking values 0 or 1) for each trial. The remaining variables represent 
the number of treatments, nt, and the number of studies, ns; r[,l] and n[,l] are the numerators and denom- 
inators for the first treatment, r[,2] and n[,2] are the numerators and denominators for the second listed 
treatment, t[,l] and t[,2] are the treatment number identifiers for the first and second listed treatments, and 
na[] is the number of arms in each trial. Text is included after the hash symbol (#) for ease of reference to the 
original data source. Note that the data are provided as a list and a table component. Both components need 
to be loaded into WinBUGS for the program to run. 

# Data (Statins example) 
list(ns=19, nt=2) 



t[,l] 


t[,2] 


na[] 


r[,l] 


n[,l] 


r[,2] 


n[,2] 


x[] 


# 


ID 


name 


1 


2 


2 


256 


2223 


182 


2221 


1 


# 


1 


4S 


1 


2 


2 


4 


125 


1 


129 


1 


# 


2 


Bestehorn 
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1 


2 


2 


0 


52 


1 


94 


1 


# 


3 


Brown 


1 


2 


2 


2 


166 


2 


165 


1 


# 


4 


CCAIT 


1 


2 


2 


77 


3301 


80 


3304 


0 


# 


5 


Downs 


1 


2 


2 


3 


1663 


33 


6582 


0 


# 


6 


EXCEL 


1 


2 


2 


8 


459 


1 


460 


1 


# 


7 


Furberg 


1 


2 


2 


3 


155 


3 


145 


1 


# 


8 


Haskell 


1 


2 


2 


0 


42 


1 


83 


1 


# 


9 


Jones 


1 


2 


2 


4 


223 


3 


224 


0 


# 


10 


KAPS 


1 


2 


2 


633 


4520 


498 


4512 


1 


# 


12 


LIPID 


1 


2 


2 


1 


124 


2 


123 


1 


# 


13 


MARS 


1 


2 


2 


11 


188 


4 


193 


1 


# 


14 


MAAS 


1 


2 


2 


5 


78 


4 


79 


1 


# 


15 


PLAC 1 


1 


2 


2 


6 


202 


4 


206 


1 


# 


16 


PLAC 2 


1 


2 


2 


3 


532 


0 


530 


0 


# 


17 


PMSGCRP 


1 


2 


2 


4 


178 


2 


187 


1 


# 


18 


Riegger 


1 


2 


2 


1 


201 


3 


203 


1 


# 


19 


Weintraub 


1 


2 


2 


135 


3293 


106 


3305 


0 


# 


20 


Wscotland 



END 

# Initial values 

# Initial values for delta can be generated by WinBUGS. 
#chain 1 

list(d=c( NA, 0), mu=c(0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0), B=0, sd=l] 
#chain 2 

list(d=c( NA, -1), mu=c(-3,-3, 3,-3,3, -3,3,-3,3,-3, -3,-3,3,3,-3, 3,-3,-3,3), B=-l, sd=3] 
#chain 3 

list(d=c( NA, 2), mu=c(-3,5,-l,-3,7, -3,-4,-3,-3,0, 5,0,-2,-5,1, -2,5,3,0), B=1.5, sd=0.5) 

Program 1(b): Rinomial Likelihood, Logit Link, Fixed Effects, Meta-Regression with Subgroups (Statins 
Example) 



# Binomial likelihood, logit link, subgroup 

# Fixed effects model with one covariate 
model! 
for(i in l:ns){ 

mu[i] ~ dnorm(0,.0001) 
for (k in l:na[i]) { 

r[i,k] ~ dbin(p[i,k],n[i,k]) 

# model for linear predictor, covariate effect relative to treat in arm 1 

logit(p[i,k]) <- mu[i] + d[t[i,k]] - d[t[i,l]] + (beta[t[i,k]]-beta[t[i,l]]) * x[i] 
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * aog(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k]))) 



# *** PROGRAM STARTS 

# LOOP THROUGH STUDIES 

# vague priors for all trial baselines 

# LOOP THROUGH ARMS 

# binomial likelihood 



resdev[i] <- sum(dev[i,l:na[i]]) 



# summed residual deviance contribution for this trial 



totresdev <- sum(resdev[]) 
d[l] <- 0 
beta[l] <- 0 
for (k in 2:nt)i 

d[k] ~ dnorm(0,.0001) 



# Total Residual Deviance 

# treatment effect is zero for reference treatment 

# covariate effect is zero for reference treatment 

# LOOP THROUGH TREATMENTS 

# vague priors for treatment effects 
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beta[k] <- B # common covariate effect 

} 

B ~ dnorm(0,.0001) # vague prior for covariate effect 

} #*** PROGRAM ENDS 

# Initial values 
#chain 1 

list(d=c( NA, 0), mu=c(0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0), B=0] 
#chain 2 

list(d=c( NA, -1), mu=c(-3,-3,3,-3,3, -3,3,-3,3,-3, -3,-3,3,3,-3, 3,-3,-3,3), B=-l) 
#chain 3 

list(d=c( NA, 2), mu=c(-3,5,-l,-3,7, -3,-4,-3,-3,0, 5,0,-2,-5,1, -2,5,3,0), B=1.5) 



The results (including the model fit statistics^^) of the 2 separate analyses and the single analysis using 
the interaction model for fixed and random treatment effects models are shown in Table A3. For the FE 
models, convergence was achieved after 10,000 burn-in iterations for separate analyses (20,000 iterations 
for the joint analysis), and results are based on 50,000 samples from 3 independent chains. For the RE 
models, 40,000 burn-in iterations were used for the separate analyses, 50,000 burn-in iterations were 
used for the joint analysis, and results are based on 100,000 samples from 3 independent chains. Note 
that in an FE context, the 2 analyses deliver the same results for the treatment effects in the 2 groups, 
whereas in the RE analysis, due to the shared variance, treatment effects are not quite the same: they are 
more precise in the single analysis, particularly for the primary prevention subgroup, where there was 
less evidence available to inform the variance parameter, leading to very wide Cris for all estimates in 
the separate RE meta-analysis. However, within the Bayesian framework, only the joint analysis offers a 
direct test of the interaction term p, which, in both cases, has a 95% CrI that includes the possibility of 
no interaction, although the point estimate is negative, suggesting that statins might be more effective in 
secondary prevention patients. 



Table A3. Statins Example^^: Posterior Mean, SD, and 95% Credible Interval (CrI) for the Log-Odds Ratio 
(LOR) and Odds Ratio (OR) for All-Cause Mortality When Using Statins for Primary and Secondary 

Prevention Groups 

Fixed Effects Random Effects 

Primary Prevention Secondary Prevention Primary Prevention Secondary Prevention 



Separate Analyses Separate Analyses 





Mean 


SD 


CrI 


Mean 


SD 


CrI 


Mean/ 
Median 


SD 


CrI 


Mean/ 
Median 


SD 


CrI 


LOR 


-0.11 


0.10 


(-0.30, 


-0.31 


0.05 


(-0.42, 


-0.18 


0.74 


(-2.01, 


-0.36 


0.16 


(-0.72, 








0.09) 






-0.21) 






1.12) 






-0.06) 


OR 


0.90 


0.09 


(0.74, 


0.73 


0.04 


(0.66, 


1.12 


3.65 


(0.13, 


0.71 


0.11 


(0.49, 








1.09) 






0.81) 






3.07) 






0.94) 


a 


NA 


NA 


NA 


NA 


NA 


NA 


0.79 


0.98 


(0.06, 


0.16 


0.23 


(0.01, 




















3.90) 






0.86) 


resdev 


16.9" 






29.0'> 






11.9" 






28.3'' 






Pd 


6.0 






15.0 






9.3 






16.8 






DIG 


22.9 






44.0 






21.1 






45.1 







(continued] 
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Table A3. (Countinued) 





Single Analysis with Interaction Term, 


P, for Subgroup 


Single Analysis with Interaction Term, 


P, for Subgroup 
















Mean/ 






Mean/ 








Mean 


SD 


CrI 


Mean 


SD 


CrI 


Median 


SD 


CrI 


Median 


SD 


CrI 


p 


-0.21 


0.11 


(-0.42, 








-0.29 


0.26 


(-0.86, 














0.01) 












0.20) 








LOR 


-0.11 


0.10 


(-0.30, 


-0.31 


0.05 


(-0.42, 


-0.07 


0.20 


(-0.48, 


-0.36 


0.16 


(-0.72, 








0.09) 






-0.21) 






0.36) 






-0.07) 


OR 


0.90 


0.09 


(0.74, 


0.73 


0.04 


(0.66, 


0.95 


0.21 


(0.62, 


0.70 


0.11 


(0.49, 








1.09) 






0.81) 






1.43) 






0.94) 


o 


NA 


NA 


NA 


NA 


NA 


NA 


0.19 


0.20 


(0.01, 


























0.76) 








resdev'^ 


45.9 












42.6 












pD 


21.0 












24.2 












DIG 


66.9 












66.8 













Posterior median, standard deviation [SD), and 95% credible interval (CrI) of the between-trial heterogeneity (a) for All-Cause Mortality When Using 
Statins for the fixed and random effects models and measures of model fit: posterior mean of the residual deviance (resdev), number of parameters (pj, 
and deviance information criterion (DIC). LOR <0 and OR <1 Favor Statins. NA, not applicable. 
'^Compare w^ith 10 data points. 
•^Compare with 28 data points. 
■^Compare with 38 data points. 



Example 2. Certolizumab: Continuous Covariate 

This example and results are described in the main text. The WinBUGS code for RE meta-regression model 
with a continuous covariate and noninformative priors is given in program 2(a), and the FE code is given in 
program 2(b). The relevant nodes to monitor are the same as in program 1. 

Program 2(a): Rinomial Likelihood, LogitLink, Random Effects, Meta-Regression witli a Continuous Covariate 
(Certolizumab Continuous Covariate Example) 

# Binomial likelihood, logit link, continuous covariate 

# Random effects model for multi-arm trials 

modeli #*** PROGRAM STARTS 

for(i in 1 :ns)i # LOOP THROUGH STUDIES 

w[i,l] <- 0 # adjustment for multi-arm trials is zero for control arm 

delta[i,l] <- 0 # treatment effect is zero for control arm 

mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines 

for (k in l:na[i]] ( # LOOP THROUGH ARMS 

r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood 

# model for linear predictor, covariate effect relative to treat in arm 1 (centring) 

logit(p[i,k]) <- mu[i] + delta[i,k] + (beta[t[i,k]]-beta[t[i,l]]) * (x[i]-mx) 
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k]]]] 

} 

resdev[i] <- sum(dev[i,l:na[i]]) # summed residual deviance contribution for this trial 

for (k in 2:na[i]] { # LOOP THROUGH ARMS 

delta[i,k] ~ dnorm(md[i,k],taud[i,k]) # trial-specific LOR distributions 

md[i,k] <- d[t[i,k]] - d[t[i,l]] + sw[i,k] # mean of LOR distributions (with multi-arm trial correction) 
taud[i,k] <- tau *2*(k-l)/k # precision of LOR distributions (with multi-arm trial correction) 

w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,l]]) # adjustment for multi-arm RCTs 
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sw[i,k] <- sum(w[i,l:k-l])/(k-l) # cumulative adjustment for multi-arm trials 



totresdev <- sum(resdev[]) 
d[l]<-0 
beta[l] <- 0 
for (k in 2mt)( 

d[k] ~ dnorm(0,.0001) 

beta[k] <- B 



# Total Residual Deviance 

# treatment effect is zero for reference treatment 

# covariate effect is zero for reference treatment 

# LOOP THROUGH TREATMENTS 

# vague priors for treatment effects 

# common covariate effect 



B ~ dnorm(0,.0001) 

sd ~ dunif(0,5) 
tau <- pow(sd,-2) 



# vague prior for covariate effect 

# vague prior for between-trial SD 

# between-trial precision = (1/between-trial variance) 

# *** PROGRAM ENDS 



The data structure is the same as in program 1, but now we have more than 2 treatments being compared 
and need to add the mean covariate value mx to the list data, for centering. 
# Data (Gertolizumab example - covariate is disease duration) 
list(ns=12, nt=7, mx=8.21) 



t[,l] 


t[,2] 


na[] 


n[,l] 


n[,2] 


r[,l] 


r[,2] 


x[] 


# 


ID 


Study name 


1 


3 


2 


63 


65 


9 


28 


6.85 


# 


1 


Kim 2007 (37) 


1 


3 


2 


200 


207 


19 


81 


10.95 


# 


2 


DE019 Trial (36) 


1 


3 


2 


62 


67 


5 


37 


11.65 


# 


3 


ARMADA Trial (34) 


1 


2 


2 


199 


393 


15 


146 


6.15 


# 


4 


RAPID 1 Trial (40) 


1 


2 


2 


127 


246 


4 


80 


5.85 


# 


5 


RAPID 2 Trial (41) 


1 


5 


2 


363 


360 


33 


110 


8.10 


# 


6 


START Study (57) 


1 


5 


2 


110 


165 


22 


61 


7.85 


# 


7 


ATTEST Trial (51) 


1 


5 


2 


47 


49 


0 


15 


8.30 


# 


8 


Abe 2006 (50) 


1 


4 


2 


30 


59 


1 


23 


13.00 


# 


9 


Weinblatt 1999 (49) 


1 


6 


2 


40 


40 


5 


5 


11.25 


# 


11 


Strand 2006 (62) 


1 


7 


2 


49 


50 


14 


26 


0.92 


# 


12 


GHARISMA Study (64) 


1 


7 


2 


204 


205 


22 


90 


7.65 


# 


13 


OPTION Trial (67) 



END 

To estimate all the pairwise log-odds ratios, odds ratios, and absolute effects, for covariate values 0, 3, and 
5, vector z could be added to the list data as list[z=c[0,3,5), nz=3). 

# Initial values 

# Initial values for delta can be generated by WinBUGS. 
#chain 1 

list(d=c( NA, 0,0,0,0,0,0), mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0), sd=l, B=0) 
#chain 2 

list(d=c( NA, -1,1,-1,1,-1,1), mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3), sd=0.5, B=-l) 
#chain 3 

list(d=c( NA, 2,-2,2,-2,2,-2), mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0, 5, 0), sd=3, B=5) 



In a RE model with Uniform(0, 5) prior for a, the heterogeneity parameter is not identifiable. This is 
because there is a trial with a zero cell and not many replicates of each comparison. Due to the paucity of 
information from which the between-trial variation can be estimated, in the absence of an informative prior 
on a, the relative treatment effect for this trial will tend toward infinity. We have therefore used an informa- 
tive half-normal prior, represented by the solid line in Figure Al, which ensures stable computation: a ~ 
Half-Normal(0,0.32"). 
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d ~i 1 1 1 I ~ i " r~ 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 



This prior distribution was chosen to ensure that, a priori, 95% of the trial-specific ORs he within a factor 
of 2 from the median OR for each comparison. Under this prior, the mean a is 0.26. To fit the RE meta- 
regression model with this prior distribution, the line of code annotated as "vague prior for between-trial SD" 
in program 2(a) should be replaced with the following 2 lines: 

sd ~ dnorm(0,prec)l(0,) # prior for between-trial SD 
prec <- pow(0. 32,-2) 

This prior should not be used unthinkingly. Informative prior distributions allowing wider or narrower 
ranges of values can be used by changing the value of prec in the code above. 

In this example, the posterior distribution obtained for a is given by the dotted line in Figure Al and shows 
that the range of plausible values for a has not changed much, but the probability that a will have values 
close to zero has decreased. 

Program 2(b): Rinomial Likelihood, Logit Link, Fixed Effects, Meta-Regression with a Continuous Covariate 
(Certolizumab Continuous Covariate Example) 

# Binomial likelihood, logit link 

# Fixed effects model with continuous covariate 

modeli # * * * PROGRAM STARTS 

for(i in 1 :ns)i # LOOP THROUGH STUDIES 

mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines 

for (k in l:na[i]) { # LOOP THROUGH ARMS 

r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood 

# model for linear predictor, covariate effect relative to treat in arm 1 

logit(p[i,k]) <- mu[i] + d[t[i,k]] - d[t[i,l]] + (beta[t[i,k]]-beta[t[i,l]]) * (x[i]-mx) 
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k]))) 

} 

resdev[i] <- sum(dev[i,l:na[i]]) # summed residual deviance contribution for this trial 

} 
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totresdev <- sum(resdev[]) # Total Residual Deviance 

d[l] <- 0 # treatment effect is zero for reference treatment 

beta[l] <- 0 # covariate effect is zero for reference treatment 

for (k in 2:nt){ # LOOP THROUGH TREATMENTS 

d[k] ~ dnorm(0,.0001) # vague priors for treatment effects 

beta[k] <- B # common covariate effect 

} 

B ~ dnorm(0,.0001) # vague prior for covariate effect 

} #*** PROGRAM ENDS 

# Initial values 
#chain 1 

list(d=c( NA, 0,0,0,0,0,0), mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), B=0) 
#chain 2 

list(d=c( NA, -1,1,-1,1,-1,1), mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3), B=-2) 
#chain 3 

list(d=c( NA, 2,-2,2,-2,2,-2), mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0, 5, 0), B=5) 



Example 3. Certolizumab: Baseline Risk 

This example and results are described in the main text. The WinBUGS code for the meta-regression model 
with adjustment for baseline risk for random and fixed treatment effects is similar to programs 2(a) and 2(b), 
respectively, but now x[i] is replaced with mu[i] in the definitions of the linear predictor. The variability of 
the normal prior distribution needs to be reduced to avoid numerical errors (this only minimally affects the 
posterior results). The relevant nodes to monitor are the same as in program 1. 

Program 3(a): Binomial Likelihood, Logit Link, Random Effects, Meta-Regression with Adjustment for 
Baseline Risk (Certolizumab Baseline Risk Example) 

# Binomial likelihood, logit link 

# Random effects model for multi-arm trials 

modeli #*** PROGRAM STARTS 

for(i in l:ns)i # LOOP THROUGH STUDIES 

w[i,l] <- 0 # adjustment for multi-arm trials is zero for control arm 

delta[i,l] <- 0 # treatment effect is zero for control arm 

mu[i] ~ dnorm(0,.001) # vague priors for all trial baselines 

for (k in l:na[i]) { # LOOP THROUGH ARMS 

r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood 

# model for linear predictor, covariate effect relative to treat in arm 1 

logit(p[i,k]) <- mu[i] + delta[i,k] + (beta[t[i,k]]-beta[t[i,l]]) * (mu[i]-mx) 
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k]))) 

} 

resdev[i] <- sum(dev[i,l:na[i]]) # summed residual deviance contribution for this trial 

for (k in 2:na[i]) { # LOOP THROUGH ARMS 

delta[i,k] ~ dnorm(md[i,k],taud[i,k]) # trial-specific LOR distributions 

md[i,k] <- d[t[i,k]] - d[t[i,l]] + sw[i,k] # mean of LOR distributions (with multi-arm trial correction) 
taud[i,k] <- tau *2*(k-l)/k # precision of LOR distributions (with multi-arm trial correc- 



tion) 



w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,l]]) # adjustment for multi-arm RGTs 
sw[i,k] <- sum(w[i,l:k-l])/(k-l) # cumulative adjustment for multi-arm trials 
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totresdev <- sum(resdev[]) 
d[l]<-0 
beta[l] <- 0 
for (k in 2:nt){ 

d[k] ~ dnorm(0,.0001) 

beta[k] <- B 



# Total Residual Deviance 

# treatment effect is zero for reference treatment 

# covariate effect is zero for reference treatment 

# LOOP THROUGH TREATMENTS 

# vague priors for treatment effects 

# common covariate effect 



B ~ dnorm(0,.0001) 

sd ~ dunif(0,5) 
tau <- pow(sd,-2) 



# vague prior for covariate effect 

# vague prior for between-trial SD 

# between-trial precision = (1/between-trial variance) 

# *** PROGRAM ENDS 



The data structure is the same as example 1, but without variable x[]. 
# Data (Gertolizumab, baseline risk) 
list(ns=12, nt=7, mx=-2.421) 



t[,l] 


t[,2] 


na[] 


n[,l] 


n[,2] 


r[,l] 


r[,2] 


# 


ID 


Study name 


1 


3 


2 


63 


65 


9 


28 


# 


1 


Kim 2007 (37) 


1 


3 


2 


200 


207 


19 


81 


# 


2 


DE019 Trial (36) 


1 


3 


2 


62 


67 


5 


37 


# 


3 


ARMADA Trial (34) 


1 


2 


2 


199 


393 


15 


146 


# 


4 


RAPID 1 Trial (40) 


1 


2 


2 


127 


246 


4 


80 


# 


5 


RAPID 2 Trial (41) 


1 


5 


2 


363 


360 


33 


110 


# 


6 


START Study (57) 


1 


5 


2 


110 


165 


22 


61 


# 


7 


ATTEST Trial (51) 


1 


5 


2 


47 


49 


0 


15 


# 


8 


Abe 2006 (50) 


1 


4 


2 


30 


59 


1 


23 


# 


9 


Weinblatt 1999 (49) 


1 


6 


2 


40 


40 


5 


5 


# 


11 


Strand 2006 (62) 


1 


7 


2 


49 


50 


14 


26 


# 


12 


CHARISMA Study (64) 


1 


7 


2 


204 


205 


22 


90 


# 


13 


OPTION Trial (67) 



END 

# Initial values 

# Initial values for delta and other variables can be generated by WinBUGS. 
#chain 1 

hst(d=c( NA, 0,0,0,0,0,0), mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0), sd=l, B=0) 
#chain 2 

hst(d=c( NA, -1,1,-1,1,-1,1), mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3), sd=0.5, B=-l) 
#chain 3 

Hst(d=c( NA, 2,-2,2,-2,2,-2), mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0, 5, 0), sd=3, B=5) 



Program 3(b): Rmomial Likelihood, Logit Link, Fixed Effects, Meta-Regression with Adjustment for Rasehne 
Risk (Certohzumab Rasehne Risk Example) 

# Binomial likelihood, logit link 

# Fixed effects model with one covariate (independent covariate effects) 
modeli #*** PROGRAM STARTS 
for(i in 1 :ns){ # LOOP THROUGH STUDIES 

mu[i] ~ dnorm(0,.001) # vague priors for all trial baselines 

for (k in l:na[i]) { # LOOP THROUGH ARMS 

r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood 

# model for linear predictor, covariate effect relative to treat in arm 1 

logit(p[i,k]) <- mu[i] + d[t[i,k]] - d[t[i,l]] + (beta[t[i,k]]-beta[t[i,l]]) * (mu[i]-mx) 
rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators 
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dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution 
+ (n[i,k]-r[i,k]) * aog(n[i,k]-r[i,k]] - log(n[i,k]-rhat[i,k]]]] 



resdev[i] <- sum(dev[i,l:na[i]]] 



# summed residual deviance contribution for this trial 



totresdev <- sum(resdev[]) 
d[l] <- 0 
beta[l] <- 0 
for (k in 2:nt){ 

d[k] ~ dnorm(0,.0001) 

beta[k] <- B 



# Total Residual Deviance 

# treatment effect is zero for reference treatment 

# covariate effect is zero for reference treatment 

# vague priors for treatment effects 

# common covariate effect 



B ~ dnorm(0,.0001) 



# vague prior for covariate effect 

# *** PROGRAM ENDS 



# Initial values 
#chain 1 

list(d=c( NA, 0,0,0,0,0,0), mu=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0), B=0) 
#chain 2 

list(d=c( NA, -1,1,-1,1,-1,1), mu=c(-3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3), B=-2) 
#chain 3 

list(d=c( NA, 2,-2,2,-2,2,-2), mu=c(-3, 5, -1, -3, 7, -3, -4, -3, -3, 0, 5, 0), B=5) 
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